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NUMERICAL SOLUTION FOR CHEMICALLY GENERATED 


WAVES IN A DILUTE, ISOTHERMAL ATMOSPHERE 
I. INTRODUCTION 

A study is made of the behavior of chemically generated waves in a simplified 
atmosphere. The atmosphere is assumed unbounded, isothermal, one-space- 
dimensional and initially quiescent. At an initial time a dissociation reaction, 
AB + J- A + B + J, commences and drives the subsequent wave motion. The 
fraction of reactant in the atmosphere, , is assumed to be small. The system 
of governing equations is then expanded in terms of the small parameter, X^j, 
and an asymptotic integral solution as X^ 0 is obtained. 

The anal 3 i;ical development is presented in Part 1 (Eberstein and Shere, 1971) 

» 

and is briefly summarized below. 

The non-dimensionalized system of equations was expanded about the parameter 
Xq , so that each dependent variable is represented in a series of the form 

00 

f (t,z) = 2 ^ f<N)(t,z)XN (1.1) 

N= 0 

where higher order terms may be dropped as Xq -* 0; t is a nondimensional time 
and z is a nondimensional altitude. Dropping terms of second order and above, 
one obtains: 

f (t, z) = f(0) + (1.2) 

Further define subscripted quantities for the density p, fraction of reactant 
dissociated a and temperature T such that 


1 


/ 


f(«> = f, 


f(l) f(0) f 


and for the velocity u such that 


=U„s. 


It follows that 


f = f(0) (1 +Xo f.,0, 


In terms of subscripted quantities, the linearized, non-dimensionalized system 
of equations describing the atmosphere is given below; 


^^( 1 ) ^^( 1 ) 
B t ^ d z 


- U-, . = 0 


1 ^^(1) , ^^(1) T 

gBt 3z 3z 


^^( 1 ) , 1 ^'^( 1 ) . .. _ 

3t 7- 1 Bt ^ Bt 


= 1 - exp [-kp t] 


( 1 . 10 ) 


The above equations are respectively, continuity, momentum, energy, and extent 
of reaction. 

S is defined as; 


B - (B/To - 1) 


( 1 . 11 ) 






m 


^ o' 


I'UiV.' 


where B is the enthalpy of reaction; 

Tq is unperturbed temperature; 
pQ is unperturbed density; 
g is gravitational acceleration. 

Defining 6 = e "‘'/2 and taking the following initial conditions for 

6>(0, z) = 0 

6^ (0, z) = - (7 - 1) Bkp exp (- 3 z/2) 

(0, z) = (7 - 1) Bk^ exp (- 5 z/2) 

One obtains the solution 0 = 6 ^ + 9 ^ where 

(t, z) = [- (7 _ 1) Bkp/(c yT)] (t exp(-c t /J)] exp (-3 z/2) 
and 9 ^ (t, z) satisfies the equation: 

— — M [0] = f (a) exp (- z/2) 

Bt2 


with M [ . ] the operator 




Integrating (1.14) with respect to time twice, one obtains 

M[0] =w(t,z) - M[0j] 

where 

w (t, z) = - kp (7 - 1) B - 


( 1 . 12 ) 


(1.13) 


(1.14) 


(1.15) 
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e+*. 


n, = exp [- (5z/2 + k_ t e"*)3 


Then, Q (t, z) is given as follows: 


t,z)= r 

•'O 


(t, z) = W (t, z, t) dr. 

h 


Taking ^ = tc, the expression forW becomes; 


W(t,Z,T)= J Jo /(tc - j 


Q (t, z, tj) d 7j 


where 


Q(r, z, 7j) = [w (t, z +7?) + w (t, z -tj)]/2c^ 


The first order terms are then: 


T (t, z) = 2 )\ 


u(t, z) = 


- f T 
~ Jz. 


Ba 


dz 


where z. is at the earth's surface. 


Defining 


E =. 


7-1 


Ba 


(1.16) 

(1.17) 

(1.18) 

(1.19) 

( 1 . 20 ) 
( 1 . 21 ) 

( 1 . 22 ) 
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and substituting into (1.21) yields: 

Bu _ BE 
B z B t 

and 


(1.23) 



-E(t.z) 


since E (0, z) = 0. 

From the linearized solution we have: 


/0(t, z) = 



(1.24) 


From (1.21) - (1.23), 


ft 

udt 

Jo 



(1.25) 


Reversing the order of integration one obtains: 





E(t, z)dz. 


(1.26) 


substituting (1.26) into (1.24) yields the expression: 


p(t, z) = E(t, z) - 



E(t, z) dz . 


(1.27) 
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II. COMPUTING 


The mathematical development of the solution was briefly outlined above. A 
discussion of numerical procedures is now presented. 

The core of the numerical problem is evaluation of the integral expressions 


02 (t.z) = 



W (t, z, r) dr 


( 2 . 1 ) 


and 


W (t, z, r) = 




Q(t, z, 17) dr? 


( 2 . 2 ) 


Evaluation of the above integrals was done numerically using Euler's method: 


N 

0(t,z)= ^ Wj(t,z,rj)Ar 

i :: 0 


(2.3) 


where Ar = t/N andr. = iAr. 
Similarly, 






z, Vj)^v 


(2.4) 


where 


A?7 = (t - t) c/M 
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and 


7j. = j At? 

The Bessel function was evaluated using an integral expression. Accuracy of the 
routine was checked against tables in the handbook of Chemistry and Physics. 
Agreement was five significant figures or better. 

The size of the integration mesh was decreased until the results became invari- 
ant to further decreases of mesh size. 

The above procedure works quite well, except for situations where the time 
becomes very long. If the integration steps are kept small, then computer time 
becomes long. Conversely, if the number of intervals is kept constant, then , 
accuracy suffers. Also, one would like to observe how the perturbation profile 
develops with time. It would thus be desirable to break up the integral for 02 
into a series of integrals as follows: 


K .Atk 

(t,z)= > W(t,,z,r) 

kTi ‘'At(k-l) 


dr 


( 2 . 6 ) 


where K is chosen such that 

t = KAt 

The above procedure could be followed quite arbitrarily if W were not a function 
of final time, t. At this point it becomes helpful to invoke some physical reasoning. 


Let us first examine the expression for W, i.e. 
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W(t,Z,T) = 


Q(t, z, T7)d7? 


Jo - 7?2'j 

The quantity \/ (tc - rc)^ - 17^ is familiar from the theory of wave propagation, 
(t - t)c is the distance that a wavelet has travelled in time (t - r). Thus the 17 
integration goes to c (t - t), or to the limits of the physical region being consid- 
ered, provided that there is no reflection at the boundaries. 

Q (t, z, 77) = [w (r, z + rj) + w (t, z - 77) ]/ 2 c^ 

where w = 0 outside the reaction zone. Let emax be the distance from the 
furthest point in the region considered to the furthest point in the reaction zone. 
If tc > emax, one has summed all the contributions. Defining At > emax/c, it 
is thus permissible to write: 


tc -Tc 




e 


2 


pAt - 2At -nAt 

Wdr + I Wdr + • • • I Wdr 
Jo JAi J(ti-l)At 


with 


W = 


ic(kAt-T) 

Jo 


and 


(k-l)At^r<;k At. 

The quantity t - r ranges from zero to At making the values of the Bessel 
function independent of k. Since the range of 77 is thus tc, ic becomes possible 
to compute a matrix of Bessel functions 
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MJo(I. J) 


r. 

where the index I refers to (t - t ) (I) and the index J refers to 77 (J). Q, however, 
depends on the actual time elapsed, and must thus be computed for each value of 
k. The integral for W thus becomes 


where 


W(K. I) = 


JMAX 

z: 


j= 1 


MJo(I.J)Q(K.J) Atj(J) 


r = r(K) 


Also, 


I MAX 

^ 2 = 

I 1 


and 


(I) = 


KMAX 



W (K, I) Ar (K) 


The rate of change of temperature, or 0, is needed to evaluate the velocity. This 
is computed using the forward difference approximation; 


dl ~ " t - t 

n m 
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The derivatives of 9 ^ and a are evaluated from exact analytical expressions. 
Once the temperature and alpha derivatives have Leen evaluated, the velocity is 
estimated by a simple numerical integration in Z , taking the velocity at the 
earth's surface to be zero. The accuracy of this integration can be improved by 
taking a finer mesh of Z. Once the velocity is known it becomes a simple matter 
to evaluate the density and pressure deviations. 

It should be noted that the step sizes for Z were between 0.1 H and 0.5 H. Thus 
the velocity and density estimates are generally less reliable than the tempera- 
ture estimates. 

III. NUMERICAL RESULTS 

At very short times a pulse is seen propagating up and down from the reaction 
zone. The initial pulse has the appearance of a discontinuity. Eventually the 
pulse passes outside the range of the computation regime, and a pseudo-steady 
pattern is established in which the qualitative behavior of the parameters does 
not change. However, the quantitative values increase to a maximum, then 
decay. The development and decay of the pseudo-steady patterns may be seen 
in Figures 1 through 3. 

The type of behavior observed may be partially explained by analogy with a shock 
tube whose driven end is semi-infinite. Initially the shock passes, then a pseudo- 
steady state is established and eventually decays. The above analogy is incom- 
plete, since, the atmosphere also behaves like an elastic medium resulting in the 
establishment of something like a standing wave pattern. However, an acoustic 
treatment would be incorrect. Firstly, the gravity restoring force is important. 
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and secondly, the reaction generates pulses all of one sign, either compression 
or rarefaction. Perhaps another analogy is a spring with weights at the end. 

When weights are suddenly added or removed, wave patterns are set up in the 
spring. Now, consider that a series of weights are added or removed in succes- 
sion. Also, let the spring be very stiff near the bottom, becoming progressively 
more elastic toward the top. 

But enough analogy. Let us proceed to examine the patterns in Figures 1-4. 

Figure 1 shows buildup of the 6 profile. The reaction zone was chosen to be 
0.5 H deep, the reaction was exothermic, and the nondimensional rate constant 
was 0.01. The reference value, was 10"^ sec~^ ; thus time is given in units 
of 100 seconds. The characteristic time of the reaction is defined as the point 
where kpt = 1. Strictly speaking, one needs kpp^°^ t = 1, but = 1 at the 
bottom of the reaction zone. Since drops off exponentially with altitude, 
while kp remains constant, it follows that the characteristic time increases ex- 
ponentially with altitude until the end of the reaction zone is reached. It may be 
shown that 

kpt = kp* t* 

where the starred quantities are dimensional. For kp = 10"^ sec"^ the charac- 
teristic time is 10^ in non-dimensional units, or 10^ sec = 2.8 hrs. At z = 0.5 H, 
the characteristic time becomes 4.6 hrs. 

The maximum value of 0 is at z = 0 and grows with time. A second maximum 
is found at z = 10, this secondary maximum also grows with time. Essentially, 

9 follows a bessel function type of altitude pattern as might be expected. The 
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rates of growth of 0 at z = 0 and at z = 10 are shown in Figures 2 and 3. The 
growth is seen to be almost linear initially, tapering off to a maximum, then de- 

ca 5 dng. Figure 4 shows development of the temperature profiles. 

1 

Outside the reaction zone the velocity equation becomes: 

^ _ - 1 ^ 

Bz y - 1 Bt 

During the buildup of temperature we have 

BT 

everywhere except in a small region near six scale heights. One may conclude 
that velocity is generally negative during the buildup phase, especially above 
8 scale heights. This conclusion is indeed correct. 

Physically we know that subsidence in an isothermal atmosphere results in 
heating. Conversely, an upward motion gives rise to cooling. This is precisely 
the kind of behavior observed, so we may conclude that the result is physically 
consistent. Similarly, convergence of velocity results in compression while 
divergence gives rarefaction. These combined effects are illustrated in Figure 5 
which, incidentally is for an endothermic reaction, five scale heights deep. 

If we had an adiabatic lapse rate, then physical argument would lead us to expect 
disappearance of the most pronounced part of the temperature wave. However, 
the velocity wave would not disappear; consequently, a temperature effect would 
become visible at an altitude where the lapse rate became less than adiabatic. If 
the atmosphere has an inversion, then the temperature wave would be especially 
pronoimced. 
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Several special cases will now be considered: 


The first case to be considered is a severe rain storm. Take saturated air at 
5 km altitude (270 degrees Kelvin) and condense 25% of it to ice, using a reaction 
characteristic time of 3 hours and a reaction depth of half of a scale height. The 
mixing ratio, is 0:002, and the non-dimensional enthalpy, B, is a little less 
than 200. The reaction is of course exothermic. Since the extent of a rainstorm, 
or even a hurricane, is not large enough to really warrant a one-dimensional 
approximation, the storm was assumed to be approximated by a line source 
centered at 6 H below ground level. The above approximation results in a 4 
amplitude fall-off relative to the pure one dimensional case. It is realized that 
the two-dimensionality correction employed is quite arbitrary, and that a genuine 
, two-dimensional solution of the atmospheric equations of motion should really be 
used. However, the results, shown in Figure 4, seem to agree reasonably well 
with experimental observations, as shown by Eberstein and Theon (1971). 

Our calculations were taken to ten scale heights above the bottom of the reaction 
zone. Whereas ten scale heights was chosen quite arbitrarily, there are nonethe- 
less compelling reasons for limiting the vertical extent to which computations 
are carried out. Firstly, the one-dimensional assumption becomes ever less 
meaningful as the vertical extent of space is increased. Secondly, the tempera- 
ture perturbation involves an exponential in altitude, i.e., 

1=6 e ^/2 

with the consequence that small errors in 6 can give rise to large temperature 
errors as Z becomes large. Also, the non-dissipation and isothermal atmosphere 
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assumptions lose validity as one considers effects propagating over large 
distances. A more detailed and comprehensive theory would be needed to study 
the effect of severe thunderstorms or hurricanes on regions in the ionosphere 
and above. But we can say at this point is that such effects would definitely be 
expected. Such a conclusion is indeed borne out by experimental observations. 
Thus Bauer (1958) has shown a convincing correlation between hurricane passage 
and electron concentration in the Fj layer of the ionosphere. More recently, 
Davies and Jones (1971) have reported association between ionospheric distur- 
bances in the region and severe thunderstorms. Davies and Jones believe 
that the ionosphere is perturbed by infrasonic disturbances generated by mechani- 
cal motions of the thunderstorms. However, the above authors do not believe that 
these disturbances are due to buoyancy oscillations. We would suggest that heat 
released by the storm induces vertical motion (rising or subsiding, but not rapidly 
oscillating in time) thus influencing the electron concentration and transmission 
properties of the Fj layer. 

While the ionosphere is outside our quantative reach, the ozonosphere is relatively 
accessible, and will now be discussed. 

Reed (1950) suggests qualitative explanations in terms of vertical and horizontal 
motions for correlations between ozone concentration and weather phenomena. 

The existence of an ozone-weather relationship is described as well known. Reed 
specifically considers subsidence at high altitudes as one of the means by which 
ozone concentration is increased. 

Our model predicts that a severe storm will cause considerable subsidence at 


ozone altitudes. 


i 


At the time of computation we had considered the velocity and density information 
to be of secondary importance. A rather crude Z mesh was thus used to save 
computer time, with the result that the quantative velocity profiles must be con- 
sidered approximate. Nonetheless, we are confident that a severe storm causes 
large and sustained subsidence at ozone altitudes. It would be very interesting 
to use a fine vertical mesh computation to obtain a quantitative estimate of the 
ozone concentration change. 

One might also consider what effect changes in the ozone layer might have on the 
rest of the atmosphere. A reaction having the thermal properties of ozone dis- 
sociating to molecular oxygen was considered. The reaction characteristic time 
was taken to be 20 minutes. This results in a maximum perturbation at approxi- 
mately 4 hours. The type of perturbation profile attained is shown in Figure 5. 

A 2 degree Kelvin cooling at the first maximum seems quite reasonable (Krueger, 
1971). It is then found that at the second maximum (9 H above the bottom of the 
reaction zone) the temperature change is 15 degrees. The associated density 
perturbation is 6%, and vertical velocity is 40 cm/sec at 10 minutes, going to 
5 cm/sec at 4 hours. Considering that the second maximum is above 90 km 
where large atmosphere variation are frequently found, one must conclude that 
upper atmosphere effects of ozone variations are not very important. 

A very different conclusion is reached when it comes to the upper atmosphere 
effects of aurorae. Insofar as a rather large amount of heat is rapidly released 
in a small altitude regime, the thermodynamic effect of aurorae is very similar 
to that of a severe thunderstorm or hurricane. The type of behavior expected is 
thus generally similar to that shown in Figure 4. If one takes an initial heating 
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rate of 25 ergs/cm^-seo and a reaction characteristic time of 20 minutes, with 
heating concentrated in half a scale height, then the quantitative deviations are 
approximately the same as those shown in Figure 4. If the mean heating rate is 
25 ergs/cm^-sec, then the deviations become twice as large. According to 
D, Heath (1971) auroral heating rates vary between 10 ergs/cm^-sec and 100 
ergs/cm^-sec with heating concentrated in less than one scale height. Corre- 
sponding characteristic times vary between 10 minutes and 100 minutes. Since 
our solution is linear, it becomes possible to estimate upper atmosphere effects 
anywhere in this range. The maximum temperature deviation predicted is then 
in the order of 500 degrees at some 200 km for the case of a mean heating rate 
of 100 ergs/cm^-sec. The actual value of 500 degrees must, of course, be taken 
with a large grain of salt. However, the theory does predict a large temperature 
increase well above the main auroral display altitude. 

IV . CONCLUSIONS 

A one-dimensional model for impulsive heat release in the atmosphere has been 
developed. The theory described is intended as a simple tool to study the effects 
of impulsive heat release. Such heating, or cooling, is found to cause large dis- 
turbances at higher altitudes. 

The next important development would be to include a second space dimension. 

It is recognized that solutions including more than one space dimension do exist. 
Best known among these are acoustic waves, gravity waves, and tidal waves. 
However, all three above mentioned waves involve special restricted solutions 
to the atmospheric equations. Specifically, acoustic and gravity waves have 
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sinusoidal space and time behavior. Tidal waves have a sinusoidal time behavior 
and a spatial behavior described in terms of Hough functions. The above theories 
are quite good for evaluating the long distance propagation of periodic distur- 
bances. However, these theories may not readily be employed to study the short 
distance effects of impulsive heat releases in the atmosphere. Our one- 
dimensional theory has been an initial step toward an analytical solution to the 
problem of impulsive heat releases in the atmosphere. Thunderstorms, hurri- 
canes, chemical reactions, and aurorae have been discussed as important natural 
sources of impulsive heat release. 


17 


REFERENCES 


Bauer, S. J., "Correlation between Atmospheric and Ionospheric Parameters", 
Geofisica Pura e Applicata, Vol. 40 (1958/11) pp. 235-240. 

Bauer, S. J., "An Apparent Ionospheric Response to the Passage of Hurricanes", 
Journ. Geophysical Research, Vol. 63 (1958) pp. 265-269. 

Davies, K. and Jones, J. E., "Ionospheric Disturbances in the Fj Region Asso- 
ciated with Severe Thunderstorms", J. Atmos. Sci., Vol. 28 (1971) pp. 254-262. 

Eberstein, I. J. and Shere, K., "A Linearized Approach to Chemically Generated 
Waves in a Dilute, Isothermal Atmosphere", Goddard X-document X-651-71- 
71, February 1971. 

Eberstein, I. J. and Theon, J., "Estimated Effect of the Release of Latent Heat 
in the Troposphere on Vertical Temperature Structure in the Atmosphere" 
submitted to Journal of Atmospheric Sciences, 1971. 

"Handbook of Chemistry and Physics", 43rd Edition, 1961-62. Chemical 
Rubber Publishing Co., Cleveland, Ohio. 

Heath, D. Personal Communication, Goddard Space Flight Center, Jan. 1971. 

Krueger, A. Personal Communication, Goddard Space Flight Center, April 1971. 

Reed, R. J., "The Role of Vertical Motions in Ozone-Weather Relationships", 

J. of Meteorology, Vol. 7, pp. 263-267, (1950). 


18 







J 



ALTITUDE: SCALE HEIGHTS 


1 

j 

■i 

i 

I 

, ] 

!t 

i 

TEMPERATURE DEVIATION FOR STORM . I 

3hr HALF-LIFE AT Z = 1h ! 


1 Rnn QPP 








